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We calculate the local density of states (LDOS) for an infinite graphene sheet with a single 
centro-symmetric out-of-plane deformation, in order to investigate measurable strain signatures on 
graphene. We focus on the regime of small deformations and show that the strain-induced pseudo- 
magnetic field induces an imbalance of the LDOS between the two triangular graphene sublattices 
in the region of the deformation. Real space imaging reveals a characteristic six-fold symmetry 
pattern where the sublattice symmetry is broken within each fold, consistent with experimental and 
tight-binding observations. The open geometry we study allows us to make use of the usual contin¬ 
uum model of graphene and to obtain results independent of boundary conditions. We provide an 
analytic perturbative expression for the contrast between the LDOS of each sub-lattice, showing a 
scaling law as a function of the amplitude and width of the deformation. We confirm our results by 
a numerically exact iterative scattering matrix method. 

PACS numbers: 72.80.Vp,73.23.-b,72.10.Fk,77.80.bn 


Introduction .— Graphene under strain has been largely 
discussed in the literature and explored for different ge¬ 
ometries, with particular features providing alternative 
routes to confine and control its charge carrier^®. Sig- 
nificant development in the theoretical description of 
strained graphene elucidated how its electronic proper¬ 
ties are modified on strained surfaces. At the microscopic 
level, a general deformation is described by modifications 
in the atomic positions which reflects in the Hamilto¬ 
nian as local changes in the hopping parameteJ 4 3. In 
the continuum model these changes appear as an effec¬ 
tive gauge field, and electrons with momentum around 
the Dirac valleys move in the deformed region as in the 
presence of a pseudomagnetic fielcP. Strain also pro¬ 
duces a deformation potential, i.e., a scalar field simi¬ 
lar to a local chemical potential that can affect electron 
dynamics 6 ! Very recently, measurements in high-quality 
graphene samples on particular substrates suggested a 
strong connection between random fluctuations in strain 
and transport propertied. 

The use of strain effects to engineer graphene elec¬ 
tronic properties has also been explored in several ex¬ 
periments in the last yeard^. As one of the most 
relevant findings, Levy et al. were able to show the 
presence of pseudo Landau levels generated by giant 
pseudomagnetic fields induced by homogeneous strain in 
graphene nanobubbled. This experimental confirmation 
that strain can have striking effects on the electronic 
properties of graphene has been followed by other ex¬ 
periments that explore the effect of different geometries 
as a path to control graphene electromechanical!}! 10 13 . 
A generic deformation of a graphene sheet can cause in¬ 
homogeneous strain, which results in an effective non- 
uniform pseudomagnetic field and provides an experi¬ 


mental test-bed to explore the interplay between highly 
tunable magnetic fields and Dirac fermions. For exam¬ 
ple, a scanning tunneling microscope (STM) tip has been 
used not only to probe samples, but also to continuously 
deform graphene nanomembranes, demonstrating elec¬ 
tronic confinement due to non-uniform pseudomagnetic 
fieldJ 10 . For a similar experimental setup, Mashoff et al. 
obtained atomically resolved STM images of stable and 
lifted regions of graphene 11 . Whereas a hexagonal ar¬ 
rangement of the carbon atoms was found at unstrained 
regions, as expected for monolayer graphene 1 ^, within the 
strained area a triangular pattern of bright spots was ob¬ 
served, signaling a symmetry breaking between A and B 
sublattices in some regions. At the time the authors spec¬ 
ulated the effect to be caused by an instability in which 
the different sublattices acquire a zigzag configuration 
with respect to the substrate. However, a local sublat¬ 
tice rearrangement requires energies that are prohibitive 
within the regime of STM imaging making such scenario 
rather unlikely. Atomistic tight-binding models 16 18 have 
predicted the development of such asymmetry but a con¬ 
tinuum, symmetry-based description remains missing. 
Similar patterns were observed in earlier works but re¬ 
mained unexplained alsd^. These results indicate an 
incomplete understanding of the fundamental electronic 
properties of graphene samples where local manipulation 
produce effective inhomogeneous gauge and scalar fields. 

In this work, we approach this problem by investigat¬ 
ing the electronic properties of a graphene sheet in the 
presence of an axially symmetric out-of-plane deforma¬ 
tion. The strain produced by such distortion is repre¬ 
sented by a pseudomagnetic field with trigonal symmetry 
and embodies a good approximation to standard exper¬ 
imental configurations, while still allowing for analytical 


2 


treatment. We use a scattering formalism based on the 
continuum description of graphene to address the ques¬ 
tion of confinement of electrons due to this deformation. 
In particular we calculate the LDOS and show that a no¬ 
ticeable imbalance in the distribution of charge density 
between the two graphene nonequivalent sublattices ap¬ 
pears even for small deformations, providing a possible 
explanation for the experimentally observed sublattice 
asymmetry. We perform exact numerical calculations 
and show that these results are well described within an 
analytical perturbative approach for small deformations. 
We analyze the dependence of the maximum LDOS con¬ 
trast between sublattices on energy and strain strength, 
providing a scaling dependence with the parameters of 
the deformation. Finally, we also show that the effective 
scalar field introduced by strain minimally modifies the 
predicted sublattice asymmetry. 

Model. — The electronic properties of undeformed 
graphene are, for low energies and large system sizes, 
governed by two copies of a two-dimensional (2D) Dirac 
Hamiltonian Ho = upp • cr where up ~ 10 6 m/s is the 
velocity of graphene electrons, p the electronic momen¬ 
tum around the K (K’) point, and cr = (a x , a y ) are Pauli 
matrices reflecting the pseudospin degree of freedom as¬ 
sociated with the sublattice structure of the honeycomb 
lattice 14 . The strain is produced by a mechanical defor¬ 
mation modeled with a height-profile h( r) that is centro- 
symmetric and is written generically as h( r) = A ho(r/b), 
where ho contains the radial profile, and the parameters 
A and b describe amplitude and effective radius of the de¬ 
formation. In the following, to illustrate our results, we 
consider the case of a Gaussian bump with height pro¬ 
file ho(x) = e~ x . Note however that our results below 
are qualitatively valid for a generic profile ho with axial 
symmetry. 

The effect of such deformation on the electronic prop¬ 
erties in the continuum limit are described within the 
theory of elasticity. For an out-of-plane deformation the 
strain tensor of elasticity^ is derived from the height 
profile h according to e^- = \d{hdjh, which in polar co¬ 
ordinates (r, 0) reads 



FIG. 1: (Color online) Schematic view of the graphene lattice 
with a magnified out-of-plane deformation, (b) Pseudomag- 
netic field created by a deformation with a Gaussian height 
profile as in (a), (c) Spatial profile of the LDOS for sublat¬ 
tice A in the presence of a bump, see Eq. (13). Bright (dark) 
spots indicate an increase (decrease) of LDOS compared to 
the undeformed graphene sheet. For sublattice B, the effect 
is exactly opposite. 


For the radial symmetric deformation, the associated 
pseudomagnetic field B = V x A shares the trigonal 
symmetry of the graphene lattice, 

B z {r) = ab 0 (r/b) sin 39 (3) 

where the spatial profile is given by the function bo(x) = 
— f'(x). For the Gaussian-shaped deformation, one 
has bo(x) = Sx 3 e~ 2x2 . 

In addition to the gauge field, the electrons are also 
exposed to a scalar potential proportional to the trace 
of the strain tensor. In our model it is represented by 
V(r) = —g s af(r/b ) with g s = 3eV as a typical valued 
The low-energy electronic properties in the presence of 
the deformation are hence described by 


e = af(r/b) 


cos 2 0 sin 0 cos 0 \ 
sin 6 cos 0 sin 2 9 J 7 


(1) 


where a = A 2 /b 2 sets the strength of the strain, while its 
spatial distribution is contained in the function f{x) = 
\ [h'o(x )] 2 . For the Gaussian profile, one has f{pc) = 
2x 2 e~ 2x2 . 

In the presence of the deformation, electrons experi¬ 
ence the strain as a gauge field 

A(r >—^“ /<r/6) ( “in*)' (2) 

where we chose the zigzag direction to lie along the x- 
axiiP. The coupling constant is g v = fShv^/ZcL « 7eV^, 
being (3 = |<91og£/<91oga| ~ 3, and t and a the hopping 
parameter and the lattice constant of graphene. 


H = v f [p + eA(r)] • cr + V(r) (4) 

In this article we consider a bump that is smooth on 
the scale of the lattice constant, such that a coupling 
between the valleys can be neglected 16 18 20 21 . Moreover, 
we consider an infinite graphene sheet containing a single 
deformation, hence our results are independent of finite 
size effects and boundary conditions^ 

Perturbation theory .— In this section we present ana¬ 
lytic results for the change in the LDOS produced by 
the scattering of electrons off the deformation obtained 
with a perturbative approach in real space. We consider 
therefore small deformations that allow for an expansion 
in the parameter a. 

From now on we set h = v-p = 1, and work around 
the K valley. The effect of the K’ valley is discussed at 
the end of this section. We split the Hamiltonian in the 
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kinetic part Hq and the perturbation V, 

V(r) = eA(r). < r=^ + ° (r) A f r ^, (5) 

where we defined A±(r) = e(A x (r) db iAy{r)) = 

(^Aa^) a f{ r /b)e T2te - We neglect the scalar potential in 
this part, we will include its effect in the next section. 

Let us start with the states of the Dirac equation in 
the absence of the deformation. Here, we take circular 
waves as a set of basis states, 


are interested in the different sublattice occupations, we 
further introduce the sublattice-resolved LDOS 

va / b (£, r) = ^(4' m (r)|7 :, A/B|^m(r)) (11) 


where V A /b are projectors on the respective sublattice 
A/B. For undeformed graphene, evaluating Eq. (11) with 
the free states |<L^(r)) produces the well-known value of 


,(°) 


[ J l—i/2|( CT )] 2 = 


5 

47T ’ 


( 12 ) 


I^M) 


I ^ ^imO f e l 6 , 2 J\ m -i/2\{er) \ 

V 4 tt \ isgn(m)e ie / 2 J| m+1/2 |(er) )' 


(6) 

where 6 denotes the energy of the Dirac fermions (which 
we assume to be positive, for simplicity), and m is a 
half-integer index labeling the states according to their 
angular momentum. J n (x) denotes the Bessel function 
of n-th order. We chose a normalization such that 



Our goal is to find the scattering state |T m (r)), that 
replaces |<l>^(r)) when the bump is present. This is de¬ 
termined by the Lippmann-Schwinger equation 

l^(r)> = |*£>(r)> + J dr'G(r, r / )V(r , )|\t , TO (r / )) (8) 

which contains the Green’s function of graphene, 


G(r,r') = -in 

rri 


|^ ) (r))($y ) (r')|, r < r' 

|$& ) (r))($<° ) (r')|, r>r 


(9) 


where we defined 


l*£°(r)> 



isgn(m)e ie / 2 H^ +1/2l (er) ) 


( 10 ) 

= J M ( x ) =b iY^{x) are Hankel functions of first 
and second kind. A derivation of Eq. © is given in the 
Supplementary Material 22 . 

For small deformations, we solve the Lippmann- 
Schwinger equation in the Born approximation, and re¬ 
place the scattering state |T m (r)) on the right-hand side 


of the equation by the unperturbed state |<L^(r)). Note 
that our perturbative approach is valid for g s , v a e. 
The explicit form of the resulting scattering states is 
shown in the Supplementary Material 22 . The trigonal 
symmetry of the pseudomagnetic field underlies the cou¬ 
pling between angular momentum states differing by 3. 

The LDOS is obtained by calculating v(e, r) = 
^) m (T m ( r ) |T m (r)). The new states are properly nor¬ 
malized to leading order in <a, as the linear in a correc¬ 
tion is orthogonal to the unperturbed state. Since we 


for the LDOS per sublattice. 

We now want to discuss the effect of the deformation 
on the LDOS. Specifically we address the limit eb <C 1, 
which is the relevant case for experiments with a radius 
of a few lattice constants. In this case, one may simplify 
the results by using the asymptotic expressions of the 
Bessel and Hankel functions for small arguments r < 
Upon retaining only the leading contribution for small 
energies, one finds for the corrections 5va,b = ^b-^ab 
to leading order in a 

Su A (e,r) 5v B (s,r) / 3A 2 

v A (e,r) vb(£,v) ba 

with the function 


sin 3 6g(r/b) (13) 


9{x) = ~ 3 / dy y 3 f (y) (14) 

x Jo 


To leading order in <a, one can replace by v a ,b in 


the denominator of Eq. (13). Notice that the relative 


LDOS correction has no dependence on energy. Thus, 
the deformation changes the local occupation in the dif¬ 
ferent sublattices in opposite directions , and their spatial 
distribution shares the symmetry of the underlying pseu¬ 
domagnetic field with a radial distribution governed by 
the function g(x). Specifically for a Gaussian height pro¬ 
file, one finds 


»<*> - i 


-2nd 


(1 T 2x 2 T 2x^) 


(15) 


The spatial distribution of the change in LDOS for sub¬ 
lattice A according to Eq. (13) is shown in Fig. [l] (c) for 
a Gaussian deformation of typical dimensions. 

A quantity of experimental relevance is the LDOS con¬ 
trast between sublattices, defined as 


C = 2 


|Fa ~ vb\ 

Z'A + "B 


(16) 


which is plotted as a function of the radial distance in 
Fig. ! (a). Note that, for a fixed width b of the de¬ 


formation, Eq (13) indicates that the contrast C scales 


quadratically with the amplitude of a centro-symmetric 
deformation. This scaling is shown for a Gaussian defor¬ 
mation in Fig. [3] and compared with the exact numerical 
results presented in the next section. 
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2 r/b 3 


a real magnetic field that has the same sign in both val¬ 
leys). These two effects ensure that their contributions 
to the sublattice occupancy contrast are identical. 
Numerics .— In this section we discuss briefly our ex¬ 
act numerical approach for the continuum model, which 
is not restricted to the case of small amplitude de¬ 
formations. The results obtained confirm our findings 
described in the previous section in the corresponding 
regimes. We use a slight modification of the method in¬ 
troduced in Ref. |24j which allows for the calculation of 
the scattering matrix S for an arbitrarily-shaped scalar 
potential. The extension to include a vector potential is 
straightforward. To calculate the LDOS integrated over 
a certain (arbitrarily chosen) volume V per sublattice, 
we include a fictitious additional scattering potential 


( £ n 

o N 

lx/ 1 ' 

reV 

V 0 

£B j 

1 lo, 

else 


FIG. 2: (Color online) LDOS contrast C as a function of 
distance from the bump’s center for fixed angle 0 = 7t/2 
(A = 0.1 nm,b = 0.5 nm). (a) Comparison of C between per¬ 
turbative (solid line, blue), and exact numerical approaches 
(red points) for £ = 0.5eV. (b) Different data sets obtained 
numerically for £ = 0.5eV, g s — 0 (blue), e = 0.5eV, g s — 3eV 
(red), e = 0.1 eV,g s = 3eV (green), and e = 0.9 eV,g s = 3eV 
(yellow). 


Such potential locally changes the electron’s energy in 
V by a magnitude £a/b in the different sublattices A/B. 
The LDOS per sublattice integrated over V is then found 
from the scattering matrix S via 




s 'ir 


(18) 


£ a/b —'0 



A[A] 


FIG. 3: (Color online) Scaling of contrast C as a function 
of amplitude for fixed Gaussian width b = 0.5 nm. Curves 
obtained with perturbative (solid line, blue), and exact nu¬ 
merical methods (points, red) for e = 0.5eV are compared. 


To conclude this section let’s discuss the role of the two 
valleys K and K’ in these results. As mentioned above, 
the deformation is smooth enough that does not couple 
the valleys and their contributions add directly. To see 
that these are identical, note that the Dirac Hamiltonian 
takes the same form in both valleys when the spinors 
are written in the valley symmetric representation 23 : 
{jpAi ^Pb) around valley K, and (^b, —^a) around valley 
K’. Note that the components referring to A and B sub¬ 
lattices are interchanged between different valleys. On 
the other hand, the pseudomagnetic field enters the Dirac 
equation with opposite sign for each valley (in contrast to 


where j = A, B specifies the sublattice. Fig. [2] (a) shows a 
comparison between analytical and numerical results for 
realistic parameters. Note that in the region of small 
amplitudes the contrast obtained with the expression 
from perturbation theory follows closely the exact solu¬ 
tion given by the numerical approach. Numerical calcula¬ 
tions in the presence of the scalar potential V (r) induced 
by the deformation were carried out for different values 
of the phenomenological parameter g s . Our results, as 
shown in Fig. 2] (b), confirm that its effect on the con¬ 
trast is negligible (note that, to leading order in a, E(r) 
affects the occupation of both sublattices in the same 
way, thus not affecting the contrast (p~6])) . Furthermore, 
the energy independent value for the contrast predicted 
with the analytical approach (as long as the requirements 
e ag s , v and sb 1 are met) is also verified in the nu¬ 
merical results as shown in Fig. 1(b)- 
Conclusions. — We have shown that a centro-symmetric 
deformation, with its local breaking of the lattice trans¬ 
lational symmetry, produces a local sublattice symmetry 
breaking on the electronic properties of a graphene sheet, 
and a consequent LDOS contrast between sublattices. 
Analytic expressions within the Born approximation pre¬ 
dict the intensity of the LDOS contrast to be determined 
by the amplitude of the deformation and to be energy in¬ 
dependent for the range of validity of the approximation. 
Exact numerics carried out with scattering matrix meth¬ 
ods confirm the validity of these results, for experimen¬ 
tally realistic parameters. While our numerical approach 
allows us in principle to treat any size of deformations, 
we concentrated here on the study of small deformations 
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and showed that there is a measurable LDOS contrast 
between sublattices even in the absence of Landau lev- 
ebP. The crossover to this regime will be published else¬ 
where. Our findings here provide an alternative inter¬ 
pretation for recent experimental observations on STM 
graphene nanomembranes 25 and provide a quantitative 
way to guide the use of strain in the design of electronic 
properties of graphene flakes. 
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Appendix A: Greens functions 

In this appendix, we show a derivation of the Green 
function for graphene given in Eq. 9 in the main text. 
The defining equation reads 

(s + icr ■ V)G(r, r') = S(r — r'). (Al) 

For the derivation, it is customary to introduce the aux¬ 
iliary Green function G s ( r, r') via 

G(r, r') = (e-icr ■ V)G S (r, r') (A2) 

The advantage of the Green function G s is that it satisfies 
a scalar equation [i.e. the pseudospin degree of freedom 
has been eliminated). 

(e 2 + V 2 )G s (r, r') @= 6(r - r') (A3) 

Moreover, G s is the solution to the Helmholtz equation 
in two dimensions. We would like to express the Green 
function in polar coordinates, thus we express the Lapla- 
cian as 

A = d 2 + -dr + dg , (A4) 

and we write for the ^-function 

5(r — r') = -5(r — r') — e m(o-o ) ? (A5) 

r 27r 

n 

where the summation is running over integer fi. We then 
expand G s as 

G s (r, 9 ; r', 9') = Y r ') ( A6 ) 

l 1 


where g M (r, r') is determined by 

(e 2 +d 2 + -d r - G) gn(r,r') = -S(r - r’). (A7) 

y r r J r 

For r ^ r', is thus satisfying a Bessel differential equa¬ 
tion. We therefore write 

9 n(r,r') = a l _ l J M (sr < )H^\sr > ) (A8) 

with r> = max(r,r'), r< = min(r, r f ). In this way, we 
respect the conditions, that 

• g^ is regular at r —>> 0, 

• g^ is behaving as an outgoing wave (oc for 

r oo (this is the proper boundary condition for 
the retarded Green function 5 e + iO), 

• g^ is continuous at r = r f . 

The remaining coefficient is determined by the jump 
condition for the derivative of 

[d r g^(r,r% Z r r ,t° 0 = 1, (A9) 

which yields = |r. Thus, the result for the scalar 
Green function reads 

G s {r,9-r',9') = L ]T e^ e ~^ J| At ,(er < ) J ff | ( + ) (er > ). 

(A10) 

The Green function for graphene is then found in a 
straightforward albeit lengthy calculation from Eq. ( |A2| ) . 


Appendix B: Scattering states 

The scattering states in the Born approximation are 
obtained by inserting Eq. 5, Eq. 6 and Eq. 9 into the 
following equation 

I’Mr)) = lM(r)> + J dr'G(r,r')V(r')\<t><£\r')). (Bl) 

After insertion of the corresponding expressions, one 
finds 
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l*m(r)>= I^W) 

- ^El $ n +) (r)) [ r ' dr ' [ ^'(4 0) (r)|V(r / )|$W(r / )) 

n Jr'<r J 

~ [ r ' dr ' [ d 6 ,, ( $ n _) ( r )l V ( r, )l $ m ) ( r/ ))- 

J r /> r J 


(B2) 


The angular integration is nonzero for n = m + 3 and n = m — 3. Calculating the matrix elements, one obtains the 
scattering states given by 

p 


iw. 


with 


>( r )> =l^m ) ( r )> + 7r ^ £a ( a m(r-)l<Am+3( r )) + 6 m( r )l^m-3( r )} +Cm(r-)l<Am+3( r )) + rf m(r-)l^m-3( r )>) ( B3 ) 

a m {r) = - Sgn(m) / dr'r'f(r'/b)J\ m+5/2 \(£r / )J\ m+1/2 \{£r / ) 

Jo 

b m (r) =sgn(m-3) / dr'r'f(r'/b)J\ m _ h/2 \ (er') J|m-i/ 2 | {er') 

Jo 

poo 

Cm(r) = - sgn(m) / dr'r'f{r'/b)H^ +5/2{ (er 1 ) J\ m +i/ 2 \(^') 

Jr 

POO 

d m {r) =sgn(m — 3) / drV/CV&^+L/ai^')J| m -i/ 2 |(er'). 

J r 

(B4) 

The leading order correction to the local density of states per sublattice up to first order on <a, when the deformation 
is present, is given by 

~~ ( £ > r ) = - 7r £ £Qie3<e E J |m-l/ 2 |( er ’) ff |mb/ 2 |( er ) a ”x( r ) 

rri 

+ Tr^eae" 3 ^ E J |m-i/2|(^)^b/2|( £r ) 6 ™(r-) 


- 7r ^“ £ae3 ^ E J |m-l/2|( £r ) J |m+5/2|(e»’)Cm(?’) 


+ 7r—£<ae E ^|m-l/2|(^)^|m-7/2|(^)^m(^) + C.C, 

m 

with va = j- for clean graphene. Using the asymptotic forms of the Bessel functions for small arguments 


Jn(%) 


1 


2 n n! 






fin (%x) , n = 0 
n>1 

7 TX n ’ 


oik' finds the dominant contribution to from the second line, for m=l /2 (+c.c.): 

SUA P b • Q Of /h\ 

- =- a smSO 0 ( 770 ). 

v A a 

The LDOS correction for sublattice B is Svb/vb = —Sva/va- _ 


(B5) 

(B 6 ) 

(B7) 

(B 8 ) 
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